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1 Introduction 

A central problem in EEG is to efficiently solve the 
inverse problem. It consists of the estimation of 
location, orientation, and strength of the current 
sources in the brain with the measurements of the 
electric potential on the surface of the head. In order 
to solve this problem, we have to first obtain the 
solution of the forward problem, i.e., the calculation 
of the electric potential due to known current 
sources in a specified conducting volume. To this 
end, a variety of efficient methods have been 
developed. The most remarkable one of the existing 
methods is directly estimating potential at the 
superficial cerebral cortical surface from EEG 
recordings on the scalp with the boundary element 
method (BEM) in a realistically shaped inhomo¬ 
geneous head model. It has made great improvement 
in spatial resolution [1], Compared to other methods, 
the cortical potential imaging approach provides a 
powerful imaging means and a much-enhanced 
spatial resolution in assessing the underlying brain 
activity without modeling the electrically activated 
brain regions by current dipoles. 

In such a method, the BEM is used to account for a 
realistically shaped inhomogeneous head that 
consists of various compartments of the head, e.g., 
brain, skull, and scalp, as they can be determined 
from Magnetic Resonance Images (MRI's). It is well 
known that the BEM is an effective approach for 
calculating the forward solution. The factors 
affecting the accuracy of the BEM in the forward 
problem were studied by many researchers. Then- 
studies demonstrated that the choice of basis 
functions, density of elements, and shape of each 
element are essential for the accuracy of BEM in the 
forward problems [3-5], Though a higher element 
density can provide a more accurate representing 
surface and a more accurate solution, element 
density is limited by physical technique. The choice 
of shapes of elements is most often confined to flat 
triangular element due to its simplicity in the 
mathematical formulation of the problem and its 
capability to approximate any complex surface. 
However, the choice of basis functions is free of the 


physical constraint. The dominant choice is constant 
basis functions due to its simplicity. However, in 
this case, the numerical errors also are very high 
unless the element density is highly increased, but 
this will result in extensive computation time and 
storage demands. A better choice is linear basis 
functions. This is because in this case, the number of 
equations can be reduced about a half of that in the 
case of constant basis functions, and the accuracy 
has not decreasing in the same rate [3], An 
interesting problem is if there exists a more suitable 
choice for basis functions choice. This paper 
addresses this problem. 

In this paper, we present a new approach by means 
of BEM based on natural element interpolation 
[2,8,9] that established a more direct and accurate 
relation between cortical potential and scalp 
potential. The natural-neighbor coordinate is 
extended so that it can be used in triangular element. 
Furthermore, the equations we established are 
directly related to nodes rather than triangles, which 
can reduce the computational cost and improve 
computational accuracy. In order to validate our 
method, we also conducted computer simulations 
based on a concentric three-shell head model. 



Fig. 1 Volume V surrounded Surface S. S = S x u S 2 



Fig.2 A three-shell volume conductor model. The 
volume between S l and S 2 is denoted as V,, and that 
between S 2 and S 3 is denoted as V 2 . 



2 Methods 
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2.1 Integral Equations 

Applying Green's second identity to an isotropic 
homogeneous volume conductor, we have [8] 
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where u is the potential, r is the vector from 
observer point to source point, r - fr| , and V is the 
volume inside the surface as shown in Fig. 1. When 
no current source exists between the surfaces, 
V 2 w = 0 . Therefore, we have 

-2nu(re S) = J^«(r)V—- — Vw(r)J• dS (2) 


where uf is the potential on jth discrete point on 
surface S k ,r = |r|. 

A k j is the yth triangular element belonging to the 
surface S k . 

Assuming {h„(r)}„ =1 is a set of basis functions and 
( r m} m=i is a set of discretization points on s x or s 2 , 
we can get 


For convenience, we modify (2) to 
1 f 1 du 




K(. r m)=8 nm - ( 6 ) 

Then the unknown potential u(r) can be expanded in 
terms of h„ (r) : 


where u is the potential on the point r and 
reSjU^, dQ. = r ^dS > #$|r| , r n is an outward 

pointing vector of unit magnitude normal to surface 
element dS. 

Our knowledge on the conductivity of the brain 
tissue is very limited, which makes us have to 
assume that the head consists of a number n of 
disjoint homogeneous regions V r ■■■,¥„ with 
boundaries s l =dv i , --,s n = dv„ • Those regions are 
volume conductors with constant conductivity 
respectively, which are usually chosen to be the 
scalp, skull, cerebrospinal fluid, gray matter and 
white matter. In this paper, we only consider a three- 
shell head model that consists of scalp, skull and 
brain, respectively, as Fig.2. So scalp v t is bounded 
by a surface 5' = j' l U5,, skull v, is bounded by a 
surface s = s, U S, , and brain ;/ is bounded by a 
surface s 3 ■ 

Since no current sources exist in volume V, (scalp) 
and Vw(r) = 0 , re S t , by applying (3) to it and by 
locating the observed point on 5j or S 2 , we can get 

h( r) = —— f ud£l + —[ ndQ .—— f --^-dS. (4) 
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By discretizing the surfaces y and s, into triangle 
elements (TE's), we get the discretization points that 
are the triangles' vertices, and from (4) we have 
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u(r) = Yj*n h n( r )- < 7 ) 

By applying (7) to (5), we get 
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Combining the left-hand side of (8) with the first 
term in the right-hand side of (8), and rewriting in 
matrix format, we get 

P n U l + P U U 2 + G X2 U 2 = 0 (9) 

where u k is the column vector consisting of 
potentials at every vertices of triangular element of 
S k and P n ,P n ,G n are coefficient matrices with 
dimensions of N\ by N ] , TVj by h\ and ,v, by ,v,, 
respectively. 
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where p)={discretepoint j contained within the elementon S k }. 


Set the observer point on S 2 , by the similar steps, 
we can get 




P 2 \U l +P 22 U 2 +G 22 U 2 = 0, (13) 

Apply (3) to the volume V 2 and locate the observe 
point on s 2 or S 3 , we can get another two equations 
P 22 U' 2 + P 23 U 3 - G' 22 U' 2 + G 13 U 3 = 0, (14) 

P }2 U 2 +P 33 U 3 +G 32 U 2 +G 33 U 3 = 0, (15) 

where u 2 is the column vector consisting of 
potentials at every vertices of triangular element of 
S 2 but just inside of v 2 . The boundary condition on 
S 2 can be expressed as follow: 

U' 2 =U 2 . (16) 

Due to the constraint of the number of sensors, we 
often set the number of the discrete point on surfaces 
S 2 and S 3 that are much greater than the number of 
sensors in order to get higher accuracy. In this case, 
the coefficient matrixes all are perhaps singular. 
Therefore, we cannot get unique solutions from (9), 
(13), (14), (15), and (16). Since this ambiguity can 
be removed by deflation [6,7], so we can derive the 
unique relationship between the potential on the 
surface y and s 2 , and that between the potential on 
the surface y, and s } . The unique relations between 
scalp potential and skull potential and that between 
skull potential and cortical potential can be derived 
U 2 =AU u (17) 

U 3 =BU 2 , (18) 

From (17) and (18), we also have 

U 3 = BAU 1 = TU l . (19) 

After the decision for the basis functions is made, 
we can get the potential distribution on the cortical 
surface from the potential distribution on the scalp. 
The simple basis functions have been studied by 
several authors such as constant potential and linear 
interpolation [1,3]. 
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Fig.3 Voronoi region of an additional point 
(situation without drawn dashed) 




2.2. Natural Neighbor interpolation: Natural 
neighbor interpolation was devised by Sibson, which 
is related to triangle based interpolation [9,10], Here, 
we give an extension of it: 


We confine the interpolation in a triangular plane 
whose three vertexes are the given data point, and 
the Voronoi tessellation is the decomposition of the 
triangular plane into the cells of the nearest 
neighbourhood of the three vertexes, shown as Fig.3. 
If Sj is the Voronoi proximal cell of data point P, 
and S(x) is the proximal cell of x in a Voronoi 
tessellation of the data points plus the point P, then 
the volumes of the overlaps of S(P) with the 
proximal cells .Sj of the data points scaled by the 
volume of .S'(P) give the local coordinate of P with 
respect to the data point Pj: 


*i( P) = 


v(S(P) n S t ) 


( 20 ) 


v(A(P)) 

Therefore, the hf P) can be used to determine an 
interpolation function 

/(P) = £/ ! ,(PK. (21) 


This function and its first derivatives are continuous 
everywhere except at the data points, which makes it 
ideally suited for our method. The expressions for 
natural neighour interpolation and its first 
derivatives have been derived by Sambridge etc. 
[10], All the expressions have been used in this 
paper after a slight modification for our special case. 



Fig.4 Relative errors versus the eccentricity. The 
horizontal axis refers to the eccentricity of the 
dipole sources and the vertical axis refers to the 
relative error: NN refers to natural neighbor 
interpolation, L refers to linear interpolation, and C 
refers to constant potential method. 










3 Results 

The accuracy of the present method was tested with 
a concentric three-shell sphere head model. The 
numerical solution was compared with the values 
given by analytical expressions. In the three-shell 
sphere head model, the radii were chosen to be 1.0, 
0.92, and 0.87, respectively. The values of 
conductivity were chosen to be 1, 0.0125, and 1, 
respectively. Each surface of three-shell sphere was 
subdivided uniformly into 1642 triangles. The 
transmatrix T was calculated according (8)-(21). 
Comparisons between the methods based on natural 
neighbor interpolation, linear interpolation, and 
constant potential are made in Fig.4. The dipole 
source was located at (r-sin(^/4), 0, r cos(^/4)) 
with varying eccentricity r. 

Our results have shown that a more accurate 
estimation cortical potential was gained even when 
the linear element interpolation was used. 

4 Discussion 

The method given in this paper can provide a means 
by which a direct relationship between the SP's and 
CP's can be established. Two highly significant 
practical advantages of this method are: 1) to 
formulate the equations in term of apexes instead of 
triangles. First, if each apex corresponds to an 
electrode, we can directly determine the equations 
for the apexes from experimental data. Second, the 
number of equations can be reduced about one-half 
and the accuracy can be greatly improved especially 
when an interpolation method of high quality was 
used. 2) Natural neighbor interpolation is local and 
the derivatives of the interpolated function are 
continuous everywhere except at the data points. 
These properties meet the need for high order 
continuous basis functions. Natural neighbor 
interpolation is better than quadratic interpolation 
since the latter requires additional nodes, besides 
those at the vertices of the triangles, e.g., at the 
midpoints of the sides of the triangles. 

We did not fully investigate the extent to which the 
natural neighbor interpolation could be used to 
improve the accuracy of BEM, though the results 
have demonstrated that a high accuracy was 
achieved. The dominant shortcoming of using 
natural neighbor interpolation in BEM is that it is 
much slower in performance than linear 
interpolation. Therefore, our future work is to 
improve its algorithm implementation. 

In summary, the experiment result and theory 
analysis show that our approach is very promising to 
further our capability of exploring the relation 


between the cortical potential and the scalp potential 

noninvasively. 
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